
subroutine  compute_Hutton_Boore_ML(confidence,c,nd,stanum,evnum, &
                                    fm_stanames,nsta,nev,sigma2_HB)

    implicit none
    
    include 'seisev_EVENT_PHASE_structure.FI'
    include 'inversione_pesata_multipar.FI'
    
    integer i,k,nd,isw,nv
    integer stanum(150000),evnum(150000),nsta,nev
    integer first,last,dst,iter,n_par,lun
    real*4  dif,r,resmax,confidence,rms
    real*4  v(5000),wmag,wsd,wmax,wmin,ww(5000),ruse
    real*8  sigma2_HB,sig2_old,sig2_new,qtot,signif
    real*4  stares_vet(5000),sta_count_vet(5000)
    logical significant
    character*5 memosta
    character*5 fm_stanames(istations)
    data isw/0/

! compute all the station Hutton-Boore magnitudes    
    n_par = nsta+nev
    iter=0
    significant = .TRUE.
1   do while(significant .and. iter .lt. 100 )
!
! itero il calcolo delle magnitudo Hutton-Boore ricalcolando nuovi residui medi 
! di stazione a ogni iterazione.
! 
      iter = iter +1
    
      do i = 1,nsta
        stares_vet(i)=0.0
        sta_count_vet(i) = 0.0
      enddo


      do i=1,nd       
        r = c(i).s_h_hypdist
      ! compute station magnitude
        c(i).MLHB_sta = log10(c(i).WA_pha_amp) + 1.110 *LOG10(r /100.) +  &
                      0.00189*(r - 100.) + 3.0 + c(i).MLHB_corr
      enddo
    
!   Prima ricalcolo le magnitudo di evento con la mediana,
!   poi con la media di Huber (tanto per vedere la differenza
!   delle due varianze).
    
      first = 1
      last  = 1
      do while( evnum(last) .eq. evnum(first) )
        last = last +1
      enddo
      last = last -1
    
      sigma2_HB=0.0d00
      do while (last .le. nd  .and. first .le. nd)
    
        k= 0
        do i=first,last
          k=k+1 
          v(k) = c(i).MLHB_sta
        enddo
    
        call mediana(v,k,wmag)      
    
        k=0
        do i=first,last
          k=k+1
          c(i).MLHB_ev = wmag     !event mag
          dif = c(i).MLHB_sta - wmag
          sigma2_HB = sigma2_HB + dif*dif
        enddo   
                 
      ! vai al prossimo evento
        first = last + 1
        last = first
        do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
          last = last +1
        enddo
        last = last -1
    
      enddo        
    
      write(*,*)' calcolo delle ML Hutton-Boore con la mediana. Sigma2 = ',sigma2_HB    
    
    ! go throught the events and compute the event HB Magnitudes.
    ! ora calcolo le magnitudo HB finali con Huber

      first = 1
      last  = 1
      do while( evnum(last) .eq. evnum(first) )
        last = last +1
      enddo
      last = last -1

      sigma2_HB=0.0d00
      do while (last .le. nd  .and. first .le. nd)

        k= 0
        do i=first,last
          k=k+1 
          v(k) = c(i).MLHB_sta
        enddo

        ruse = 0.23   ! valore ottenuto da inversioni precedenti
        wmag = 0.
        wsd = 0.
        call Huber_dev(k,v,ruse,wmag,wsd,wmax,wmin,ww)
      ! memorizza i risultati   !
      !   
        k=0
        do i=first,last
          k=k+1
          c(i).MLHB_ev = wmag     !event mag
          dif = c(i).MLHB_sta - wmag
          c(i).MLHB_dif = dif
          sigma2_HB = sigma2_HB + dif*dif
          !
        ! media pesata dei residui per la stazione stanum(i):
          stares_vet(stanum(i)) = stares_vet(stanum(i)) + dif*ww(k) 
        ! somma dei pesi per quella stazione:
          sta_count_vet(stanum(i)) = sta_count_vet(stanum(i)) +ww(k) 
          
        enddo   
                 
      ! vai al prossimo evento
        first = last + 1
        last = first
        do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
          last = last +1
        enddo
        last = last -1

      enddo        


      write(*,*)' calcolo delle ML HB con Huber. Sigma2 = ',sigma2_HB    
      
      rms=SNGL(sigma2_HB)
      rms=rms/float(nd)
      rms=SQRT(rms)
      write(11,'(a,a,i3,a,f10.2,a,i5,a,f8.4,a)') &
        'Hutton-Boore',',',iter,',',sigma2_HB,',',nd,',',rms,',1.110,0.00189,3.0'
 
    ! i residui da questa iterazione secondo HB 
      do i = 1,nsta
        stares_vet(i) = stares_vet(i) / sta_count_vet(i)
      enddo
      resmax= -1.
    
      do i = 1,nd
        ! nuove correzioni (residui cambiati di segno):
        c(i).MLHB_corr = c(i).MLHB_corr - stares_vet(stanum(i))
        if(ABS(stares_vet(stanum(i))) .gt. resmax)then
          resmax = stares_vet(stanum(i))
          memosta = c(i).sta
        endif
      enddo 
      write(*,*)' residuo ABS massimo HB ',memosta,' ',resmax 
   
      if(iter.eq.1)then
        sig2_old = sigma2_HB
        significant = .TRUE.
      else
        sig2_new = sigma2_HB
        call Ftest(nd,sig2_old,n_par,sig2_new,n_par,qtot,signif)
        sig2_old = sigma2_HB
        if(signif .gt. confidence)then
          significant = .TRUE.
        else
          significant = .FALSE.
        endif
        write(*,*)' Hutton-Boore residuals; iter signif ',iter,signif
       
      endif
    enddo 
! calcolo finale della variazione dei valori di ML di stazione per ogni stazione

    do i=1,nsta
      stares_vet(i) = 0.
      sta_count_vet(i) = 0.
    enddo
    do i=1,nd
      stares_vet(stanum(i)) = stares_vet(stanum(i)) +  &
                              c(i).MLHB_dif*c(i).MLHB_dif
      sta_count_vet(stanum(i)) = sta_count_vet(stanum(i)) + 1.
    enddo
    do i=1,nsta
      stares_vet(i) = SQRT( stares_vet(i) / (sta_count_vet(i)-1.))
    enddo
    
    
! OUTPUT DEI RISULTATI HUTTON-BOORE =======================================
! *************************************************************************
! =========================================================================

    call get_lun(lun)
    open(unit=lun,file='results_HUTTON_BOORE.csv',status='unknown')
    
    write(lun,'(a)')'Hutton-Boore,iter,sigma2_HB,nd,rms,P1,P2,K'
    write(lun,'(a,i3,a,f10.2,a,i5,a,f8.4,a)') &
        'Hutton-Boore,',iter,',',sigma2_HB,',',nd,',',rms,',1.110,0.00189,3.0'

    write(lun,'(a)')'Sta-code,corr,sigma'
    do i=1,nsta
      do k=1,nd
        if(stanum(k).eq.i)then
          write(lun,'(a,2(a,f7.4))')c(k).sta,',',c(k).MLHB_corr,',',stares_vet(i)
          goto 99
        endif
      enddo
 99   continue
    enddo
    close(lun)

return
end
    
    
 